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We investigate the use of renormalisation group methods to solve partial diflFerential equations 
(PDEs) numerically. Our approach focuses on coarse-graining the underlying continuum process 
as opposed to the conventional numerical analysis method of sampling it. We calculate exactly 
the coarse-grained or 'perfect' Laplacian operator and investigate the numerical effectiveness of the 
technique on a series of 1 + 1-dimensional PDEs with varying levels of smoothness in the dynam- 
ics: the diffusion equation, the time-dependent Ginzburg-Landau equation, the Swift-Hohenberg 
equation and the damped Kuramoto-Sivashinsky equation. We find that the renormalisation group 
is superior to conventional sampling-based discretisations in representing faithfully the dynamics 
with a large grid spacing, introducing no detectable lattice artifacts as long as there is a natural 
ultra-violet cut off in the problem. We discuss limitations and open problems of this approach. 
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I. INTRODUCTION 

It is a rare event in science that a single paper con- 
tains an idea so powerful that it revolutionises an en- 
tire field. Rarer still are those gems which transform 
two or more apparently separate fields, such as the pa- 
per entitled Scaling Laws for Ising Models near Tc by 
Leo P. Kadanoff Q, published in the regrettably short- 
lived journal Physics. Indeed, this organ, edited and 
founded by P.W. Anderson and the late B.T. Matthias, 
announced as its by-line the memorable claim "An inter- 
national journal for selected articles which deserve the 
special attention of physicists in all fields" ; and per- 
haps no other journal before or since has lived up to this 
hubris. Kadanoff's famous article developed the notion 
of what came to be called "block spins" , and was cited in 
the title of the early seminal paper by K.G. Wilson which 
introduced the modern form of the renormalisation group 
RG) in both condensed matter and high energy physics 



Although these theoretical developments have become 
part of the canon of modern physics, the "spin-ofF' from 
Kadanoff's work continues to this day, as mathemati- 
cians and physicists, including Kadanoff and colleagues 
at the University of Chicago, study the singularities and 
patterns which arise in extended physical systems, gov- 
erned only by sets of partial differential equations (PDEs) 
1^. On the auspicious occasion of this 60th birthday 
Festschrift to honour Leo P. Kadanoff, it therefore seems 
appropriate to contribute a brief account of our recent 
unpublished work Q which extends the "block spin" in- 
sights and renormalisation group theory to the numeri- 
cal solution of PDEs. Some, but not all of our results 
have recently been rederived independently by Katz and 
Wiese |^ , and we shall comment on the differences in the 
sequel. This work is part of our program to utilise RG 
methods for PDEs [lljj], and is distinct from our earlier 
work applying RG iterative methods to construct simi- 



larity solutions and travelling waves 



II. THEORY OF PERFECT OPERATORS 
A. Motivation 

Discretisation is an inevitable part of numerical analy- 
sis. Let us suppose that we wish to solve a partial differ- 
ential equation numerically. The standard procedure in 
real space is to suppose that the solution u{x) is sampled 
at points Xi and an algorithm devised to approximate the 
values Ui = u{xi). If the points are equidistant with spac- 
ing dx, then we naturally require that in the continuum 
limit dx 0, the sequence Ui converges to u{x). 

The disadvantage of the sampling approach is that one 
is forced to reproduce as faithfully as possible all the 
detail and fine structure of the solution, even on a scale 
that may be of no interest or worse, beyond the regime of 
applicability of the differential equation itself. This has 
two consequences: 

• a small grid size dx must be used, which implies 
many grid points must be calculated and stored; 

• a small time step 6t is implied by the small dx, 
either for reasons of accuracy or stability of the 
numerical method. 

As a result, the numerical method is subjected to 
an unnaturally large degree of computational complex- 
ity. Finally, for a given sampling procedure there is no 
unique prescription for obtaining the equation governing 
the sampled points Ui. The only criterion for admissabil- 
ity is that it converges in the continuum limit; in practice 
of course, one seeks schemes which are numerically stable 
and attain the continuum limit rapidly. 
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From the statistical physicist's point of view, 
KadanofF's 1966 paper suggests that a more natural ap- 
proach to discarding information is to coarse-grain rather 
than to sample. Such a procedure focuses on the scale of 
interest, whilst allowing the possibility of accessing the 
continuum function to an arbitrary level of detail if de- 
sired. Suppose that we denote the coarse-graining oper- 
ator at scale A by the symbol Ca, with capital letters de- 
noting coarse-grained quantities. Then conceptually we 
need to find the operator La which connects U {X, 0) and 
U{X,t) given the microscopic time evolution operator L 
connecting 0) with u{x, t), as shown schematically in 
the commutativity diagram below: 



u{x, 0) 



Ca 



Ua{X,0) 



u{x, t) 



Ca 



UA{X,t) 



In practice, we coarse-grain onto a lattice Xi with 
a specified Ca- The freedom of choice in this coarse- 
graining operator parallels the lack of uniqueness in defin- 
ing a governing equation for sampled points Ui in the 
sampling approach; but once a coarse-graining operator 
Ca has been defined, there should be a unique prescrip- 
tion to obtain La- 

Our diagram suggests that coarse-graining and time 
evolution commute, but it is not clear that this is correct, 
even in principle. For example, an equation with a pos- 
itive Liapunov exponent might have the following prop- 
erty: two initial values ui_2(a;,0) differing only in field 
configuration at two nearby points in space, but with the 
same coarse-grained representation U{X,0), might dif- 
fer substantially at very long times t: ui{x,t) ^ U2{x,t) 
even though the coarse-grained initial fields would evolve 
identically. This example raises the interesting question 
of whether such equations are well-defined: no numerical 
procedure would be appropriate, unless the divergence of 
the trajectories was still bounded (as in a strange attrac- 
tor). 

We also need to consider the appropriate coarse- 
graining scale. Two situations are possible here. In the 
first, we suppose that the solution we wish to obtain has a 
natural scale A below which there is no significant struc- 
ture. In that case, our goal is to avoid having to over- 
discretise the problem merely in order to attain the ac- 
curacy of the continuum limit. Thus, we would like to 
be able to use as large a value for the grid spacing dx 
as possible without sacrificing accuracy. In the second 
situation, there is no such obvious scale, or at least, it is 
not known a priori, but the computational demands are 
so large that it is simply not feasible to work with a grid 
spacing dx smaller than some size A. In this case, we 
would like to minimise in some sense the artifacts that 
must inevitably arise. 

We will mainly have in mind the first situation, which 
is more straightforward because the only issue is speed of 



convergence to the continuum limit: there is no explicit 
discarding of important dynamical information. Instead 
the focus is how to remove lattice discretisation artifacts 
[^,||. We will refer to an operator or equation as be- 
ing 'perfect' if it has been constructed by coarse-graining 
appropriately so that it has no lattice artifacts; our dis- 
cussion follows the pioneering work of Hasenfratz and 
Niedermeyer in the context of lattice gauge theory ^ . 

In the second situation, however, one is making an un- 
controlled and potentially severe truncation of the correct 
dynamics, and issues of modelling must be faced. For ex- 
ample, can one model the neglected unresolved scales as 
effective renormalisations of the coefficients in the origi- 
nal PDE? Are the neglected degrees of freedom usefully 
thought of as noise for the retained large-scale degrees 
of freedom? And how can any available statistical infor- 
mation on the small-scale degrees of freedom be used to 
improve the numerical solution for the large-scale degrees 
of freedom? Such questions may perhaps be treated by 
combining constrained Monte Carlo simulation of the mi- 
croscopic degrees of freeedom with a maximum entropy 
criterion for discretisation of the large-scale degrees of 
freedom. An alternative but related approach has been 
implemented by Kast and Chorin pC[ | , who minimise the 
RMS error, estimated from knowledge of the microscopic 
probability distribution. 



B. The Perfect Laplacian 

Our goal in this section is to examine the simplest pos- 
sible problem, namely the diffusion equation in c? + 1 di- 
mensions: 



5fu(x, t) = dlu{yL,t) 



(1) 



subject to appropriate initial conditions and boundary 
conditions. 

The coarse-graining is defined with respect to a d- 
dimensional hypercubic lattice of spacing a on whose ver- 
tices reside the lattice variables ?7(n,t) defined by 



[7(n, t) = —J d'^x U{yi + na, t). 

a" J -a/2 



(2) 



The numerical scheme that we will adopt is Euler dis- 
cretisation in time, but with a perfect discretisation for 
the Laplacian: 

C/(n,i + rft) = C/(n,i) + dt AaC/(n,i) (3) 

where the perfect Laplacian at scale a is given by 

g-i^_^C/(n,t)A„;7(n,t) ^ f jj^^^i f d-'x u{x)dlu{x) 



TT(5 f C/(n,t) - 4? / ^ d'^xu{-K + na,t)\ (4) 
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where the delta function enforces the definition of the 
coarse-grained field 



U(n,t) = —: / d'^x u(:s. + na, 

a'^ J -a/2 



,t). 



(5) 



This definition of the coarse-grained operator is of course 
motivated by field theory, but is unsatisfactory for several 
reasons. First, it is admittedly ad hoc; in some sense, it is 
imposing a probability distribution on the high wavenum- 
ber degrees of freedom which is not necessarily present 
in the actual solution of the PDE being solved. Second, 
we are implicitly assuming that the coarse-graining of a 
differential equation is accomplished by simply coarse- 
graining the separate terms in the equation according to 
the prescription given above. We would prefer to able to 
start with the governing PDE and coarse-grain the entire 
equation in a systematic procedure. We hope to report 
on this methodology in a future publication. The delta- 
function constraint in Equation can be softened by 
writing 5{x) — > C(k) exp(— 2kx^), with C{k) being the 
normalisation. 

The calculation of the functional integral in Equation 
is straightforward ||^ and yields a convolution expres- 
sion 



R 



where R is a d-dimensional lattice integer, and 
P(H)^/|^e--/(p) 



n 



4sin^ Pj/2 



2Tlk 



1 

3k 



(6) 



(7) 



(8) 



Here, the vector p lies in the first Brillouin zone: \pi\ < 
TT. Evaluation of the coefficients p(R) requires numerical 
integration in general, but in the case of = 1 there is a 
simplification because 



1 4sin2p/2 ^ 



(9) 



;=-c 



{p + 2ttI)'^ ' 3k 

which can be summed by contour integration to yield 

/b)-i = icosecV2+(^-j). (10) 

In the case k = 2, 

piR) = 2Sr^o - Sba - Sb.,-1 (11) 

which looks just like the conventional lowest order finite 
difference expression for the second derivative. This is, 
however, slightly misleading, because this kernel acts on 



the coarse-grained function U{n), not the sampled value 
of the continuum function. 

The general result for arbitrary k is 



PiR) 



A 



A 



3k, 

K + 4 



1/2 



- 1 



\R\ 



if i? = 

otherwise 
(12) 



where A = —6k/ (k — 2). 

These formula are problematic to interpret for the 
cases where k < oo, because the functional integral in 
Equation (^ does not exactly enforce the definition of 
the coarse-grained field U. Thus, for example, is it consis- 
tent to coarse-grain the initial condition using Equation 
(^, whilst at the same time using the perfect Laplacian 
operator with k < oo? We will return to this issue below. 



III. NUMERICAL RESULTS IN ONE 
DIMENSION 

In this section, we present numerical results obtained 
on differential equations in one space and one time di- 
mension, but with varying character of solutions. All the 
results were initiated from the random initial conditions 
shown in Figure (0), defined on a lattice of 1024 points 
with grid spacing dx = 0.5 with periodic boundary con- 
ditions. 
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FIG. 1. Random initial condition with dx = 0.5, = 1024 
grid points. The inset shows the power spectrum of the initial 
condition, S{k). 

First we examine the diffusion equation, where high 
wavenumber behaviour dies away rapidly. As expected, 
the KG numerical method performs well compared to the 
standard algorithms. The second example is the coars- 
ening dynamics of the time dependent Ginzburg-Landau 
equation (sometimes known as model A), where the co- 
alescence of domains is the dominant behaviour; we will 
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see that over-aggressive coarse-graining of the initial con- 
ditions can lead to late time configurations which differ 
from the correct solution. Although it seems plausible 
that the statistical information is preserved, i.e. the en- 
semble averaged structure factor, we have not yet explic- 
itly checked this hypothesis. The third example is the 
Swift-Hohenberg equation, which in the parameter range 
we studied forms a lamellar phase with a well-defined pe- 
riodicity. Again, over-aggressive coarse-graining is seen 
to be counter-productive. Lastly, we studied the damped 
Kuramoto-Sivashinsky equation, which is a toy model for 
directional solidification and other systems which form 
interfacial patterns. The interesting aspect of this bench- 
mark is that a spatio-temporal chaotic phase exists for 
certain parameter ranges. As one might expect from the 
heuristic comments earlier, faithful reproduction of the 
solution is not really possible with any significant degree 
of coarse-graining. 

In order to perform these numerical experiments, we 
have had to make two uncontrolled approximations. 
First, we have followed Katz and Wiese and used the ex- 
act coarse-graining operation Equation (||) but allowed 
the Laplacian operator to have any k, not just the value 
K = oo. In particular, we have used the special value 
K = 2, where the Laplacian takes on its conventional and 
very local form. Our results are essentially unchanged 
when we use the k = oo form of the Laplacian, but there 
is a slight loss of stability of the Euler algorithm in this 
case, which sets a limit on the maximum value of dt/dx^ . 
For a general value of k, the stability limit (calculated 
for the diffusion equation) is given by 

dt 1 2 

implying that dt/dx^ < 1/6 for k = oo. 

Second, we have not yet calculated the coarse-grained 
operators corresponding to nonlinear operators such as 

and udxu: instead we have simply replaced each non- 
linear operator N{u) by writing CaN{u) = N{U). This 
approximation can be controlled in a more systematic 
derivation of the theory. 

A. Diffusion Equation 

We used the random initial condition shown in Figure 
(|l|), defined on a lattice of 1024 points with grid spacing 
dx = 0.5, and evolved it in three different ways to time 
t = 100. The first way was a benchmark calculation using 
a conventional numerical analysis Euler scheme with dt — 
0.001. 

The second way used the RG methodology we have 
described above. We coarse-grained the initial condition 
down to a smaller lattice of 128 grid points and dx = 4, 
using the coarse-graining of Equation (|^), yielding the 
function exhibited in Figure (||) . We evolved this forward 
in time using the perfect Laplacian of Equation (O) with 



K = 2, and dt = 5. Such a large value was stable because 
of the much larger value of dx than in the benchmark. 

1 I ' , ' , ' , ' , ' 1 
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FIG. 2. Random initial condition coarse-grained so that 
dx — 4 and — 128 grid points. 

In the third way, we sampled = 128 points of the 
initial condition with a uniform grid spacing of dx — 4, 
and then evolved this forward in time using the standard 
Euler method. 

The appropriate comparison is between the two calcu- 
lations using 128 points: how well does each reproduce 
the benchmark calculation? 
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FIG. 3. Comparison of results for the diffusion equation 
for t — 100 and = 128. The solid line represents the 
benchmark, the circles the coarse-grained method, and the 
triangles the conventional sampling method. 

In Figure (^) is plotted the benchmark configuration 
a.t t = 100 along with the results from the coarse grain- 
ing or sampling methods. The coarse grained result al- 
most exactly falls onto the benchmark calculation, even 
though it uses an 8-fold coarser lattice and a bigger time 
step size. Numerical stability dictates that dt < dx^ /2. 
This means that the largest time step size one can use 
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for the benchmark lattice is about dt = 0.1. Therefore, 
by coarse graining by a factor of eight and using a dt 50 
times larger than the maximum permitted for the bench- 
mark, the calculation was accelerated by a factor of 400, 
without introducing any noticeable lattice artifacts. The 
reduction in CPU time will be more pronounced when we 
go to higher dimensions. Furthermore, one sees clearly 
that with the same degree of discretisation, the uniform 
sampling result is inferior to the coarse-graining result 
and does not reproduce the benchmark calculation very 
well. 

The reason that coarse-graining works well in this in- 
stance is because our procedure preserved all relevant in- 
formation down to the coarse grained scale A. In the uni- 
form sampling method, however, one grid point is used 
to represent all those within a neighbourhood of size A; 
thus, this point is likely to be in the tail of the spectrum 
spanned by function values in the neighbourhood. 

What is an appropriate scale at which to coarse-grain? 
A should be set such that the time scale of interest is 
larger than the relaxation time T^eiax of details on a scale 
smaller than A. In the case of the diffusion equation, if 
the initial configuration is very smooth on the scale of 
A, we can safely coarse grain to that level; the inset to 
Figure (||) shows the power spectrum of the initial data, 
and reveals that the coarse-graining level used preserves 
the salient long- wavelength features. When there are sig- 
nificant high wavenumber modes, given that Treiax of a 
fluctuation on the order of dx is roughly proportional to 
dx^ , we can at most coarse grain to a level A ~ Vt with 
T being the time scale of interest. Because computa- 
tion is most consuming in simulations with long evolu- 
tion time and large lattices, we almost always want to 
coarse grain to some level. But if we are interested in the 
early time evolution of the configuration with significant 
high wavenumber modes, we should not coarse grain at 
all. 



B. Time-dependent Ginzburg-Landau Equation 

The time-dependent Ginzburg-Landau equation was 
studied with the discretisation 

u{n, t + dt) = u{n, t) -|- dt [e uin^ t) — u{n, t)^ + Aau{n, t)]. 

(14) 

It has a nontrivial fixed point u(n, t) = i\/e when 
e is positive. Unlike the diffusion equation, where the 
configuration continually flattens, one sees regions of 
u{n,t) — and u{n,t) = —\fi forming, separated by 
domain walls. As time progress, these regions coalesce 
and expand as shown in Figure (Q). 

In our simulation, e = 0.3, and the benchmark lattice 
was evolved with dx — 0.5, dt = 0.001. Coarse graining 
and sampling were used to evolve the system with = 2, 
dt = 0.1, and dx = 4, dt = 0.001. As shown in Figure 
(^, coarse graining always yields superior results than 



uniform sampling. Also, one notices that when coarse 
grained to dx = 4, there are some suprious offshoots 
in the configuration. This is because we coarse-grained 
too aggressively. The structure factor S{k) = |it(fc)p of 
the initial configuration is plotted in Figure (M. It has 
significant spectral weight beyond k = tt/S. Each coarse- 
graining step (defined as double the dx size) discards half 
of the Brillouin zone with larger wavenumber. So coarse- 
graining three times loses the modes with k > tt/S. This 
is not serious for diffusion where the dynamics is a trivial 
decay of large wavenumber modes, but is a poor approxi- 
mation for model A where domains grow from the initial 
configuration. 
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FIG. 4. Comparison of results for the TDGL equation with 
e = 0.3 for t = 100 and = 128. The solid line represents 
the benchmark, the circles the coarse-grained method, and 
the triangles the conventional sampling method. 

In Figure (||) we have compared the results when a 
smaller dx is used: with this resolution, the coarse- 
graining algorithm still captures most of the correct fea- 
tures of the solution, whereas the sampling method in- 
troduces a spurious domain. 
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FIG. 5. Same as Figure (4) but with N = 256. 
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C. Swift-Hohenberg Equation 



The discretised form of the Swift-Hohenberg equation 


u{n, t + dt) = u{n, t)+dt [e u{n, t) — u{n, t)^ 

~il + Aafu{n,t)] (15) 

was studied using the comparison methodology described 
above, and the resuhs are shown in Figure (0). 
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FIG. 7. Solution of the Swift-Hohenberg equation with 
e = 0.5, dx = 7r/2, dt = 0.001, t = 4. The solid line repre- 
sents the benchmark, the circles the coarse-grained method, 
and the triangles the conventional sampling method. 

u{n, t + dt) — u(n, t) + dt [e u{n, t) — u{n, t)\Ju{n, t) 

-~{\ + Aafu{n,t)] (16) 

where Vu(7i, t) is the usual central difference formula for 
first derivative. 
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FIG. 6. Solution of the Swift-Hohenberg equation with 
e = 0.5, dx = 7r/4, dt = 0.001, t = 4. The solid line repre- 
sents the benchmark, the circles the coarse-grained method, 
and the triangles the conventional sampling method. 

Because of the (1 -I- Aa)^ term, there is an intrinsic 
length scale of 2it and a periodic ground state. The dy- 
namics is interesting because u not only has an amplitude 
but also a phase. We choose e — 0.5. The benchmark 
lattice was evolved with dx = 7r/16, dt = 0.0001 to time 
t — A. The coarse-graining and sampling methods used 
dx = 7r/4, dt = 0.001. Again, the coarse-grained re- 
sult reproduces the benchmark quite well, better than 
sampling, but not as well as found in experiments on 
the diffusion and TDGL equations. We attribute this 
to the existence of the intrinsic length scale. When we 
coarse-grain to a level close to this period, many relevant 
modes are ignored. Indeed, when we coarse-grain further 
to dx = 7r/2, the dynamics longer follows the benchmark 
very closely, as shown in Figure (M). 



D. Damped Kuramoto-Sivcishinsky Equation 

We also considered the damped Kuramoto-Sivashinsky 
equation , which can undergo a phase transition from 
a lamellar phase to spatio-temporal chaos when the con- 
trol parameter e exceeds a critical value Cc ~ 0.68. The 
discretised equation of motion studied was 
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FIG. 8. Solution of the damped Kuramoto-Sivashinsky 
equation with e = 0.9, dx = 7r/8, dt = 0.0001, t = 16. 
The solid line represents the benchmark, the circles the 
coarse-grained method, and the triangles the conventional 
sampling method. 

We examined both the lamellar phase (e = 0.4) and 
the chaotic phase (e — 0.9). The former situation is very 
similar to that found in the Swift-Hohenberg equation. In 
the chaotic phase shown in Figure (||) , one finds that even 
coarse-graining once does not give us a result satisfacto- 
rily close to the benchmark. This is not surprising given 
that the chaotic phase has positive Liapunov exponent 
as discussed before. In Figure (H), the benchmark was 
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obtained using dx = tt/16, dt = 0.0001 to time t = 16, 
whereas the coarse graining and sampUng methods used 
dx — 7r/8, dt — 0.0002. Of course, the interesting ques- 
tion to address is whether the statistical properties of the 
dynamics are faithfuUy reproduced by coarse-graining in 
a way that is superior to samphng methods. 



IV. DISCUSSION 

We have presented a feasible alternative to the numeri- 
cal solution of differential equations, in which the contin- 
uum limit is attained from a sequence of coarse-grained 
function values rather than a sequence of sampled func- 
tions. The numerical results that we and others [^Jl0| 
have obtained are promising and a fully systematic study 
is clearly warranted. In particular, useful acceleration of 
appropriate problems is possible without loss of accuracy 
or even the need to use adaptive methods. 

There are many issues left unresolved by our work 
(which is why it has remained unpublished until this 
timely occasion). First and foremost, what is the correct 
interpretation to place on the RG scheme when k < oo? 
It can be shown that a literal interpretation requires a 
stochastic coarse-graining transformation; in this case, 
should one (in principle) average over realisations of the 
coarse-graining noise in order to obtain the appropriate 
solution of the PDE? Is it correct to take the mean of 
this distribution, which seems to offer a justification for 
using an exact (i.e. k, — oo) coarse-graining procedure 
with a K = 2 formula for the dynamics? How can one 
properly extend the philosophy espoused here to nonlin- 
ear problems in a systematic way? 

Katz and Wiese Q implicitly addressed these issues 
by deriving the coarse-grained equations of motion from 
a postulated action functional S, which they varied with 
respect to all its arguments. In fact their procedure is 
equivalent to ours, and their action S can be derived from 
our Equation (||) by making the Hubbard-Stratonovich 
transformation 



S'{[/(n,t),u(n,i),A(n,t)} = ^ / d"^ x {u{x)dlu{x) 



E 



--aM.)} (20) 



is the action (3.5) from reference g. Such a functional 
integral formulation does not possess a small parameter 
in which to make a loop expansion based around the 
classical action 



5S 



dS 



0. 



(21) 



5u dX{n,t) 

Furthermore, only by requiring the additional constraint 
dS 



dU{n) 







(22) 



can one use the coarse-grained equation of motion for any 
K consistently with the constraint of Equation (^ (eqn. 
(3.8) of reference g). 

Despite these questions, we are optimistic that the 
spirit of the program initiated by Leo P. Kadanoff in 
critical phenomena will be equally useful in the fields of 
pattern formation and partial differential equations. 
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-2kW„ 



-x^ / {8K)+ixW„ 



(17) 



with 



I fa/2 

Wn = U{n,t) -J d'^xuhi + na,t), 



(18) 



leading to the result that 



i^_^C/(n,t)A„(7(n,t) ^ f jj^ jj ^ ^- S {U [n^t) ,u{n.t) ,\{n^t)} 



(19) 



where 
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